#include <iostream>
#include <fstream>
#include <vector>
#include <string>

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include<iostream>
#include<fstream>
#include<typeinfo>

#include "util.h"
#include "diffuse.h"

// 定义方阵的维度
#define nMAX 10  

using namespace std;


int main()
{

	char str1[64], str2[64], str3[64];

	sprintf(str1, "output_data");
	if (!access(str1, 0) == 0) mkdir(str1, 0777); // if  no "data" folder, construct it for output date

	char   buf_data[40];

	int ndim = 3;
 	
	if (ndim > nMAX) { cout << "矩阵阶数超过最大数nMAX !!!" << endl;  exit(-1); }

	// double  a[nMAX][nMAX] = { 0 };//定义一个nMAX*nMAX*的矩阵，用于存放数据
	double  **a = dmatrix(1, ndim, 1, ndim); // 或者我们可以动态开辟数组，这种方式记得在程序结束后释放内存。使用方式 free_dmatrix(a, 0, nMAX, 0, nMAX);
	
	// 怎样读取数据文件，比如data.txt中存储的是矩阵A数值
	ifstream infile;//定义读取文件流，相对于程序来说是in
	infile.open("data.txt");//打开文件
	for (int i = 1; i <= ndim; i++)//定义行循环
	{
		for (int j = 0; j <= ndim; j++)//定义列循环
		{
			infile >> a[i][j];//读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
		}
	}
	infile.close();//读取完成之后关闭文件


	ifstream infile2;
	infile2.open("data_b.txt");		// 打开文件
	for (int i = 1; i <= nMAX; i++) // 定义行循环
	{
		infile2 >> b[i]; // 读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
	}
	infile2.close(); // 读取完成之后关闭文件


//////////////////////////////////// test  /////////////////////

	// Householder vector
	int n = 10;

	double *x = dvector(1, n);
	// 然后给x赋值
	double *v = dvector(1, n);

	// beta 是指向double类型的指针变量
	double *beta; 

	// 调用Householder函数
	// beta是指针变量， 取地址，可以把beta的值带出函数体
	Householder(x, v, &beta, n);

    // 给H矩阵开辟内存	 
	double **H = dmatrix(1, n, 1, n);
	zero_matrix(H, 1, n, 1, n);

	for (int i =1; i <= n; i++)
		for (int j = 1; j <= n; j++)
		{
			if (i == j)
			{
				H[i][j] = 1.0 - beta * v[i] * v[j];
			}
			else
			{
				H[i][j] = - beta * v[i] * v[j];
			}
		}
	

    return 0;
}

// 算法3.1
void Householder(double* x, double *v, double *beta, double n)   
{
	// here n is the dimension of vector x
	// Given x in R^n, compute v in R^n such that Hx=||x||_2 e_1, where H = I - beta v v^*

	double alpha = 0.0;
	double sigma = 0.0;
	for (int i = 2; i <= n; i++)
	{
		sigma += x[i] * x[i];
	}

	vec_copy(v, x);

	if (sigma == 0)
	{
		if (x[1] < 0)
		{
			v[1] = 2.0*x[1];
			beta = 2.0 / (v[1] * v[1]);
		}
		else
		{
			v[1] = 0.0;
			beta = 0.0;
		}
	}
	else
	{
		alpha = sqrt(x[1] * x[1] + sigma);   // alpha = norm2(x)
		if (x[1] < 0)
		{
			v[1] = x[1] - alpha;
		}
		else
		{
			v[1] = -sigma / (x[1] + alpha);
		}
		beta = 2.0 / (v[1] * v[1] + sigma);
	}

	return;

}


void print_vector(double *x, int n, char* buf)
{
	char buffer_x[32];
	
	FILE *fx;

	sprintf(buffer_x, "./%s/x.m", buf);

	fx = fopen(buffer_x, "w");
	for (int i = 1; i <= n; i++)
	{
		fprintf(fx, " %16.14f ", x[i]);
		fprintf(fx, "\n");
	}
	fclose(fx);

	return;
}



void print_matrix(double** phi, int n, char* buf)
{

	char buffer_phi[32];

	FILE *fphi;

	sprintf(buffer_phi, "./%s/a.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	return;
}


void print_matrix2(double** phi, double** eta, int n, char* buf)
{

	char buffer_phi[32], buffer_eta[32];

	FILE *fphi;
	FILE *feta;


	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/L.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(feta, " %16.14f ", eta[i][j]);
		}
		fprintf(feta, "\n");
	}
	fclose(feta);

	return;
}



void print_matrix3(double** phi, double** eta, double **xx,  int n,  char* buf)
{

	char buffer_phi[32], buffer_eta[32], buffer_xx[32];

	FILE *fphi;
	FILE *feta;
	FILE *fxx;

	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/L.m", buf);
	sprintf(buffer_xx,  "./%s/U.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(feta, " %16.14f ", eta[i][j]);
		}
		fprintf(feta, "\n");
	}
	fclose(feta);


	fxx = fopen(buffer_xx, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fxx, " %16.14f ", xx[i][j]);
		}
		fprintf(fxx, "\n");
	}
	fclose(fxx);

	return;
}


